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We study the renormalized perturbation theory of the single-impurity Anderson model, particu¬ 
larly the high-order terms in the expansion of the self-energy in powers of the renormalized coupling 
U. Though the presence of counter-terms in the renormalized theory may appear to complicate the 
diagrammatics, we show how these can be seamlessly accommodated by carrying out the calcula¬ 
tion order-by-order in terms of skeleton diagrams. We describe how the diagrams pertinent to the 
renormalized self-energy and four-vertex can be automatically generated, translated into integrals 
and numerically integrated. To maximize the efficiency of our approach we introduce a generalized 
fc-particle/hole propagator, which is used to analytically simplify the resultant integrals and reduce 
the dimensionality of the integration. We present results for the self-energy and spectral density 
to fifth order in U, for various values of the model asymmetry, and compare them to a Numerical 
Renormalization Group calculation. 


I. INTRODUCTION 

The Anderson Impurity Model (AIM) 1 was introduced 
in 1961 to explain the existence of localized magnetic 
moments in metals with magnetic impurities. It stands 
today as one of the best understood impurity models, 
having been intensely studied by a number of meth¬ 
ods such as the bare perturbation theorj^ED, ]/]\[ e x- 
pansiond 7 (where N deno tes the impurity orbital degen¬ 
eracy), the Bethe Ansat an d the Numerical Renor¬ 
malization Group (NRG) G® 13 (among many others; see 
Ref?'® for a review). Each of these is particularly useful 
in different scenarios: the bare perturbation theory ex¬ 
cels at weak-coupling, 1/N expansions are asymptotically 
correct for highly degenerate models, the Bethe Ansatz 
yields exact results, but only for static quantities, and 
the NRG approach cannot efficiently deal with large de¬ 
generacies. Together, these methods paint a consistent 
picture for the physics of the model, but their abundance 
highlights the challenges posed even by relatively simple 
models of strongly correlated electrons. This in-depth 
understanding of the model has served to establish it 
as a testing ground for the development of new meth¬ 
ods in many-body physics. Though the AIM is today 
well understood in the context of a single impurity, in¬ 
terest in it has been renewed in light of the Dynamical 
Mean-Field Theory 1 ®, which allows models of impurity 
lattices to be mapped — exactly in the limit of infinite 
dimensions — onto a single-impurity problem subject to 
a self-consistency consistency condition. 

An important step in the understanding of impurity 
physics was the connection made by Nozieres between 
the strong coupling fixed point of the Kondo model at low 
temperatures and Landau’s theory of Fermi liquids 1 ^ in 
which the excited states of the model are interpreted in 
terms of quasi-particles which are weakly interacting?®, 
even in the strong correlation regime. It was quickly rec¬ 
ognized that the AIM is a Fermi Liquid in all parameter 


regimes 11 "^ thoug h this idea was not fully explored until 
very recently 17 * 18 l 

In light of the Fermi liquid interpretation, the Renor¬ 
malized Perturbation Theory (RPT)P^'" 2 ' offers a conve¬ 
nient way of analysing the low-energy behaviour of the 
model. The usual Hamiltonian H of the AIM specifies 
the model in terms of the impurity level e the hybridiza¬ 
tion broadening A and the local Coulomb interaction U. 
In the renormalized theory these are replaced by effec¬ 
tive parameters id,a, A CT , U which are used to define a 
renormalized Hamiltonian H of the same form as H and 
which, by definition, incorporates the low-energy one- 
particle interactions. The propagators in the renormal¬ 
ized theory describe the quasi-particle states, making ex¬ 
plicit the one-to-one correspondence between the single¬ 
particle excitations of the non-interacting and interact¬ 
ing systems. The quasi-particle interactions can now be 
taken into account by constructing a perturbation theory 
in the renormalized parameters and organized in powers 
of U. 

The RPT has a number of appealing featured^. i n 
the presence of a magnetic field the leading term of the 
renormalized expansion is of order U , and suffices to cal¬ 
culate the zero-temperature spin and charge susceptibil¬ 
ities exactly in all parameter regimes. In the absence of 
a magnetic field the leading term in the RPT is of order 
U 2 and leads to a simple exact expression for the sec¬ 
ond derivative of the imaginary part of the self-energy at 
zero frequency. This in turn can be used to derive an 
exact expression for the T 2 coefficient of the conductiv¬ 
ity of the symmetric model in terms of the renormalized 
parameters. 

In this paper we discuss the calculation of the renor¬ 
malized self-energy using the diagrammatic RPT and 
show how this can be implemented on a computer com¬ 
pletely automatically. Our presentation is structured as 
follows: We start by describing a simple algorithm to 
generate all relevant Feynman diagrams and retain only 
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those that contribute to the renormalized self-energy. 
Though at first the RPT seems to have a more com¬ 
plicated perturbational structure than the bare theory, 
owing to the presence of counter-terms, we show how 
these can be included into the calculation with minimal 
effort by setting up the calculation in terms of skele¬ 
ton diagrams. As an intermediate step we introduce a 
simplification algorithm which dramatically reduces the 
computational complexity of the resultant integrals by 
factorising out sub-integrations that can be computed 
analytically. Finally, we carry out the numerical integra¬ 
tions and present results to order U 5 inclusive for the 
self-energy and spectral density in the strong correlation 
regime, for different values of the asymmetry, and com¬ 
pare these to results obtained using the NRG. 


II. RENORMALIZED PERTURBATION 
THEORY 

The effective Lagrangian for the Anderson model in 
the limit of an infinitely wide conduction electron band 

ipa 

C = d a ( T ) (d T - e d ,a + iA) d a (T) 

<T— T4 

+ Un t (T)ni(r). (1) 

In the bare perturbation theory one attempts to solve 
Eq. @ by taking the U = 0 (or Hartree-Fock) state as 
the starting point for a perturbation theory in powers of 
the interaction U. Though some useful information can 
be extracted this way in the weak correlation regime, the 
value of U in physically relevant systems is usually too 
large to be handled in this manner; this is precisely the 
challenge of strongly correlated physics. In particular, 
it has long been recognized that the low-energy excita¬ 
tions primarily responsible for the interesting physics of 
impurity models cannot be described perturbatively. 

The difficulty of the bare perturbation theory can be 
ultimately traced to the unfortunate choice of the U = 0 
state as the starting point for the perturbation expan¬ 
sion. As we increase U, virtual low-energy scattering 
processes will lead to the formation of quasi-particles on 
the impurity site which interact through a renormalized 
Coulomb interaction t/, and whose relation to the orig¬ 
inal particles becomes increasingly tenuous. The RPT 
seeks to address this issue by using precisely these quasi¬ 
particle states as the starting point for the perturbation 
expansion in powers of the renormalized coupling. This 
is accomplished by writing the Lagrangian in Eq. |l]) 
as a sum of a renormalized quasi-particle Lagrangian, 
which describes the quasi-particles, a renormalized in¬ 
teraction term and a remainder term, the counter-term 


Lagrangian, i.e. C = Co + Cu + £ c t, where 
£o = ^2 da(j) (d T - e d a + d a (r), (2) 

C ct = ^ da(r)(X 2 ,adr + Xi^)d a {r) + A 3 n t (r)fi|(r), 

rr=t,i 

( 3 ) 

Cu = Un^n ( 4 ) 

In this re-arrangement the term Co forms the starting 
point for the perturbative expansion and the interaction 
term is taken to be Cu + C ct . This term gives rise 
to a renormalized self-energy T, a (u ;) and two-particle- 
reducible four-vertex F||(u;i,W 2 ;w 3 ,W 4 ). The counter¬ 
terms A = (Ai jCr , \ 2 ,a, A 3 ) are then determined by impos¬ 
ing the renormalization conditions 

S ff (0) = 0, 

<9 w E CT (w) | u=0 = 0, 

f n (0,0;0,0) = t/. (5) 

The presence of the counter-terms ensures that renor¬ 
malization effects that have already been absorbed in 
the renormalized parameters are not overcounted. The 
counter-terms are real numbers; this is due to the Fermi 
liquid property of the AIM. 

In many ways, this re-organization of the Lagrangian 
in terms of the renormalized parameters is similar to 
the corresponding practice in High Energy Physics (see 
RefP^for instance). The motivation is largely the same, 
to re-express the Lagrangian in terms of physical param¬ 
eters. Nevertheless there are conspicuous differences. For 
instance, in our system there are natural upper and lower 
energy cut-offs provided by the metal’s lattice spacing 
and sample size respectively. Therefore the technicali¬ 
ties of the regularization procedure, which can be rather 
complex for many-loop calculations, do not enter our dis¬ 
cussion at all. 

For the purposes of this article we use the NRG 
to determine the renormalized parameter£p 2 l note how¬ 
ever, that in certain cases they can be determined 
entirely within RPT without appealing to an exter¬ 
nal methocP3 27 29 . Given the renormalized parameters, 
Eq. 0 is imposed in some approximation for E^ta); the 
counter-terms thus depend on the chosen approximation 
to the self-energy, whereas the renormalized parameters 
are independent of it, and are in a one-to-one correspon¬ 
dence with the bare parameters that define the model. 
We can relate the renormalized parameters e d ,a, A CT to 
the bare ones through 

£-d,a — At (erf,fr "F ^er(O)) 5 

A ct = z a A, (6) 

where 

= l-^E ff H| w= o’ (?) 
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and similarly the relate the renormalized self-energy to 
the bare quantity through 

= Zu (£ CT (w) - w0 a ,£ o .(u;)| u , = o - £, 7 ( 0 )). (8) 

We eliminate the quantities d u T, (T (a;)| w= o and £o-(0) in 
favour of id and z„ to obtain 

id,a + £<t(w) = (£<t(w) + e d ,a) ~ (1 “ (9) 

From Eqs. ©, (§ we find the renormalized interacting 
propagator 

G ct (w) = [w - e rfiCr + iAo- - £ ct (u;)] -1 = z^G^uj) (10) 

and thus deduce via a Fourier transform and_the defini¬ 
tion of the Green’s function G ct (t) = (d ( 7 (r)d CT (0)) that 

d a (r) = zl /2 d a (r) and d a (r) = z X J’ 1 ~d rj (r). 

III. AUTOMATED RPT EXPANSIONS 

In this section we will describe the automation of the 
calculation of the self-energy. For the purposes of our 
calculation we will use the T = 0 formalism. The renor¬ 
malized Green’s function which will form the basis of the 
diagrammatic expansion is 

G' 01 M =- . . ■ (11) 

w - +*A cr sign(a;) 

The T = 0 formalism has the advantage of working di¬ 
rectly on the real axis, thus avoiding the need for an an¬ 
alytic continuation, which is often fraught with its own 
numerical difficulties. More generally, for m > 1 we in¬ 
troduce the propagator 

G l ™ ] M = CM M£H MG™ M, (12) 

where £^(w) denotes the U m (strictly) term of the 
renormalized self-energy. To simplify the discussion we 
assume a zero magnetic field, so i d , a and A CT are spin- 
independent quantities. 

Our diagrammatic approach is based on the skele¬ 
ton formalism. This has the advantage of avoiding ex¬ 
plicit reference to the counter-terms Ai and A 2 . Fur¬ 
thermore, we introduce an effective interaction constant 
U e = U + A 3 that combines the A 3 counter-term with the 
renormalized interaction constant. This allows us to ten¬ 
tatively carry out the expansion of £(w) in powers of U e , 
without having to explicitly account for A 3 . Ultimately, 
our goal is to organize the calculation in powers of U, 
rather than U e . This will be discussed at the end of this 

section, but we note for now that it requires knowledge of 

- 4 

r (T _ CT (0, 0; 0, 0) order-by-order in U e , up to U e inclusive. 
We thus aim to generate and calculate diagrams for the 
self-energy and the four-vertex. 



FIG. 1: The first-order correction to the renormalized 
self-energy. The internal line corresponds to the free 
quasi-particle propagator of Eq. (Ill and the broken 
line denotes the effective interaction vertex U e . 


A. Generating the diagrams 

Consider the first-order term given by the loop dia¬ 
gram of Fig. [l] and a tree-level counter-term. Since the 
loop is frequency independent, we see, upon imposing 
Eq. that the first-order self-energy is £^(w) = 0 . 
Since T-f^O, 0; 0,0) is by definition equal to U, we find 
that A 3 = 0(U 2 ). Note that the cancellation of the 
loop to this order implies that any diagram with a loop 
will cancel with the corresponding diagram that has a 
counter-term in place of the loop, to all orders in U. For 
similar reasons we can ignore diagrams whose external 
legs attach to the same interaction vertex, for these do 
not depend on the external frequency and will cancel with 
the counter-term. In this sense, calculations in the renor¬ 
malized theory thus involve fewer diagrams than the bare 
theory. 

To calculate the self-energy as a power series in U we 
have, in principle, to explicitly include the six interaction 
vertices - five from the possibly spin-dependent counter¬ 
terms and the sixth being U — order-by-order in the 
renormalized expansion. Fortunately, in the absence of 
a magnetic field, explicit reference to the counter-term 
vertices can be avoided by setting up the perturbation 
theory in the skeleton formalism. We define a skeleton 
diagram as a diagram that does not contain self-energy 
insertions. A skeleton diagram involving n interaction 
vertices will necessarily contain 2n — 1 internal lines. By 
replacing all of them with (T 0 ! we construct a self-energy 
diagram of order ro; by replacing any 2 n — 2 lines with 
G[°l and the remaining line with GW we obtain one of 
the diagrams that contribute to £l n+ 1 l, and so forth. An 
advantage of the skeleton formalism is that it minimizes 
the number of diagrams that have to be accounted for, 
since most of the diagrams of order n can be generated 
from a lower-order skeleton diagram. Note that for n < 6 
we do not need to consider diagrams with more than one 
insertion, since a second-order diagram with two second- 
order insertions would contribute a U 6 term. 

Rather than attempt to enumerate all relevant dia¬ 
grams manually, which is tedious and error prone, we 
used the IGRAPH librarjEl to first construct all possible 
diagrams and then retain only those relevant to our calcu¬ 
lation. Since we are only interested in relatively small or¬ 
ders n — the bottleneck is the numerical integration, not 
the diagram generation — we found it sufficient to gener- 
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V 

E 

f 

2 

1 

2 

3 

2 

9 

4 

12 

58 

5 

73 

438 


TABLE I: Number of dynamic skeleton self-energy and 
four-vertex diagrams as a function of the number of 
vertices V. 


ate diagrams by simply exhaustively considering all pos¬ 
sible combinations of connecting n vertices such that no 
internal line begins and ends on the same vertex. By re¬ 
garding each Feynman diagram as a directed graph with 
multiple edges, and calculating the graph’s edge connec¬ 
tivity, we can identify and discard one-particle reducible 
diagrams. 

In graph-theoretic terms a skeleton graph can be de¬ 
fined as a graph that does not contain a subgraph isomor¬ 
phic to a lower-order skeleton graph. Consequently, to 
identify the skeleton graphs of order n we must generate 
and store all diagrams of order 31 2,3,..., n — 2. To iden¬ 
tify the non-skeleton diagrams we then bijectively map 
each (multi-edged) Feynman graph to a simple graph 
with weights encoding the s pins o f the multiple edges 
and apply the VF2 algorithirPUMl for the subgraph iso¬ 
morphism problem to the simple graphs (this ‘edge col¬ 
oring’ technique is necessary to apply the algorithm to 
multi-edge graphs). Working order-by-order we can thus 
generate all the skeleton diagrams up to our desired or¬ 
der. 

The diagrams that contribute to the two-particle- 
reducible four-vertex can be generated in a similar way. 
Whereas in the case of the self-energy we only had to 
consider the possibility that two external legs attach to 
different interaction vertices, the diagrams of the four- 
vertex admit more complex topologies which have to be 
taken into account individually. The calculation is oth¬ 
erwise identical, except for the fact that there are only 
2n — 2 internal lines corresponding to a diagram of or¬ 
der n. Our method can similarly be extended to higher 
correlation functions, but this is beyond the scope of this 
paper. In Table [T] we report the number of relevant dia¬ 
grams to fifth order in U e . 


B. Evaluating the diagrams 

Having generated all relevant self-energy diagrams we 
proceed to impose frequency conservation at each vertex. 
An n’th order diagram for the self-energy will have 2n — 1 
internal lines. These are subject to n constraints, though 
due to global frequency conservation only n—1 are inde¬ 
pendent. We thus have 2n — 1 — (n — 1) = n independent 
frequencies, each of which corresponds to an integration 
variable. We thus arrive at a linear, underdetermined, 
system of dependent equations, cuj = fl, where c is an 
n x (2n — 1) matrix encoding the constraints in the form 


y~)(in) — E(out), lj a vector with 2n — 1 components, 
each corresponding to the frequency of an internal line, 
and LI the n-dimensional vector (—f2,0,0,..., f2), where 
H is the external frequencj®. This can now be readily 
inverted to give the internal line frequencies as a function 
of the n integration variables e and a particular solution 

flp to the systenPH 

w = /e + f2 (p) . (13) 

This can be accomplished in a number of ways; we found 
it best to determine the null space using the Lenstra- 
Lenstra-Lovasz lattice basis reduction algorithnPS, as 
this results in /-matrices similar to what one would ob¬ 
tain by manually imposing frequency conservation (i.e. 
matrices whose entries are ±1 or 0). The diagrams for 
the four-vertex are treated similarly: an n’th order dia¬ 
gram will involve n—1 integration variables, resulting in 
an nx (2n —2) c-matrix and a (2n — 2) x (n — 1) /-matrix. 

It is possible, though very inefficient, to apply the 
Feynman rules at this point and proceed with the nu¬ 
merical integration. In calculations by hand, however, 
it is common to analytically factorize out any particle- 
hole or particle-particle pair-propagators. This has the 
advantage of reducing the dimension of the numerical 
integration and consequently reducing the number of in¬ 
tegrand evaluations needed to achieve a given precision. 
We can generalize this to accommodate the more general 
scenario of k (quasi-)particle/hole lines and define 


/ OO k 

dc/n G i°w+^), (i4) 

-°° i =i 

where s* = 1 for particle lines and s* = — 1 for hole lines. 
Due to the simple form of the propagator, the integral can 
be calculated analytically; the details of the calculation 
have been relegated to the Appendix. 

Having defined the fc-particle/hole propagator we will 
now describe an algorithm to identify instances of it from 
the information encoded in the /-matrix. Our strategy 
will be to first inspect the /-matrix and identify groups 
of Green’s functions that can be combined to form a 


product of the form of Eq. (141. We accomplish this by 


traversing the /-matrix column-by-column — since each 
column corresponds to an integration variable — and ex¬ 
amining each column’s non-zero entries. Our goal is ul¬ 
timately to delete from the /-matrix the columns which 
correspond to integration variables that have been ab- 
sorbed in n^ ; s, and to delete the rows that correspond to 

~ /M 

Green’s functions that comprise n^g. To avoid modify¬ 
ing the /-matrix while we are traversing it we will main¬ 
tain a list L c of columns and a list L r of rows which are to 
be removed, both of which are empty at the start of the 
search, and we will only delete the rows and columns after 
all possible simplifications have been identified. We will 
refer to the resultant simplified matrix as the g-matrix, 
to distinguish it from the original /-matrix. 

To identify instances of Il(j fc ] we scan each column m 
of the /-matrix and identify the columns with exactly 
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k non-zero entries n \,..., Hk on the rows ri,..., r^ 37 . If 
m £ L c and rq,..., rk ^ Lr then column m indeed corre¬ 
sponds to a fc-particle/hole propagator; if not, and either 
m £ L c or at least one of rq ,... rq £ Lr, it means that 
the corresponding Green’s function has already been ab¬ 
sorbed into another propagator, and we proceed to the 
next column. With the j’th fc-particle/hole propagatoi^ 
we will associate a matrix qj, which encodes its depen¬ 
dence on the remaining integration variables, and vec¬ 
tors <jj , s j describing the spin and sign configuration 


respectively. The k frequency arguments of will 

therefore be given by the rows of qj. Furthermore, we 
construct a fc-dimensional vector $7 q that contains the 
external frequency dependence of the rows rq,..., rk, i.e. 
fi q = • ■ •, In other words, the aq that ap¬ 

pear in Eq. (14) can be obtained from the i’th component 
of 


qe' + J“J q , (15) 

where e' denotes the free variables that remain after 
all propagator simplifications. After constructing the 
q—matrix we add m to L c and all of r i,..., rq to L r 
and examine the next column. When the columns of the 
matrix have been exhausted we delete all columns in L c 
and rows in L r from the original /-matrix, to obtain its 
final form g. Furthermore, to account for the fact that 
components of the original Il <p) have been absorbed into 
the various Q q we delete all the rows in L r from f2 <p) 
to obtain its final form, which we denote Finally, 

we eliminate the columns in L c from the provisional q 
matrices. 

For each diagram we thus arrive at an expression for 
the amplitude of the form 


I'D = Pd 


<\pcA 1 n( fei ) fr( fe 2) 

^ w CTl W <T2 • ■ ■ ii cri;si ii cr2;S2 


(16) 


where Pd is a complex prefactor containing powers of 
U e ,i and 2ir. In Eq. (16) it is understood that the argu¬ 


ment of Gal 1 corresponds to the V th row of 

id = ge + 


(17) 


which is similar to Eq. (13) but involves g, the simplified 


form of the /-matrix, and that the (vector) argument of 


the j’th II^j sj is given by Eq. (15), with the correspond¬ 
ing q , and Ll qj . We remark that while the final integral of 
Eq. (16) is at least one-dimensional for all self-energy di¬ 


agrams, some four-vertex diagrams factorize completely. 

To illustrate this with an example consider the diagram 
of Fig. [2] We find that the amplitude is proportional to 


ei-ea+^G^-es + tt) 

x GL 0 l( e i)GL 0 l(e 2 )GL°l(e3), (18) 



FIG. 2: One of the two diagrams that contribute to the 
self-energy to third order in U e . 


or in matrix notation, 


T 

0 

-r 


tt 

0 

1 

-i 

, n (p) = 

n 

1 

0 

0 

0 

0 

1 

0 


0 

.0 

0 

1. 


.0. 


(19) 


We begin with the first column and identify the first 
particle-particle pair propagator by the first and third 
row entries. The propagator is 11^-,s, where <x = (ct, —ct), 
because the first and third Green’s functions have spin a 
and —a respectively, and s = (1,1) because f n = / 31 = 
1. Having ‘used’ the Green’s functions in question, we 
mask their corresponding entries in / by appending ‘1’ 
to L c and ‘1,3’ to L r so that L c = {1} and L r = {1,3}. 
We temporarily associate this propagator with the ma¬ 
trices 


qi = 


1 o -1 
1 0 0 


flqi — 


( 20 ) 


constructed from the first and third rows of / and fi 
in Eq. (19). We now examine the second column of / 
in Eq. (191 and note that the two non-zero entries are 


contained in the second and fourth rows, neither of which 

( 2 ) 

is in L r . This is thus another instance of n ; , ,, ,,, 

and after updating L c = {1,2}, and L r = {1,3,2,4}. 
We associate this propagator with a q 2 matrix equal to 
Eq. (20). We proceed to the third column and examine 


the possibility that it corresponds to a triple propagator. 
We conclude that it does not, since its non-zero elements 
appear on rows 1,2, 5, the first two of which are already 
in L r . Having exhausted the columns of /, the final step 
is to remove all columns in L c and all rows in L r from / 
to obtain g, and all column in L c from q 1; q 2 to obtain 
their final simplified form and remove the references to 
the variables of integration that have been eliminated. To 
conclude the procedure we also delete all rows in L r from 
to obtain H < ' 9 \ We thus have / = [1], H^ 9 ' 1 = [0] 

and 


qi = q 2 = 


rjq 2 — fl qi — 


and the amplitude in Eq. (18) simplifies to 

f deiGi°](ei)[ng_ (7);(lil) (e 1 + H,0) 


l 2 


( 21 ) 


( 22 ) 
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The final step in our calculation is the numerical in¬ 
tegration, which is handled with adaptive quadrature 39 
and cubature method^!. By evaluating vectors of points 
at every iteration, the calculation can be seamlessly par¬ 
allelized to take advantage of modern computer architec¬ 
tures. Note that it is possible, especially in the absence 
of a magnetic field, for two distinct diagrams to be equal 
numerically. We test for this by evaluating each diagram 
individually at a randomly chosen frequency of order A, 
and, if an equal pair is found, ensuring that only one 
of the diagrams is included in the integration, with an 
adjusted prefactor, is included in the integration. A sim¬ 
ilar situation occurs in the case of particle-hole symme¬ 
try, where id = 0 , rendering G^(w) an odd function of 
frequency and leading to the numerical cancellation of 
several diagrams. 


C. Assembling the RPT 


So far we have been working order-by-order in the ef¬ 
fective interaction U e = U+ A 3 , expressing the self-energy 
as 


EM = ^ 7 n Hf/ e n . (23) 

71 = 2 

Ultimately, we aim to obtain £(w) as a power series in 
U rather than U e . To accomplish this, we note that the 
counter-term A 3 is defined as in Eq. <©• By calculating 
the renormalized four-vertex at zero frequency we can 
thus obtain U as a power series in U e , 


where 

5 2 = 72 M, 

5 3 = ~ 2 a 2 j 2 (uj) +73M, 

64 = (5a| - 2a 3 ) 72 (u) - 3a 2 7 3 (w) + 74M, 

S3 = (—14a2 + 12a 3 a 2 — 204) 72 (w), 

+ (9«2 - 3a 3 ) 73 (w) - 4a 2 'y 4 (cj) + 75 (w). (28) 

These relations enable us to deduce the renormalized ex¬ 
pansion from the bare one and show explicitly how the 
inclusion of counter-terms results in the re-organization 
of the series. 


D. Checks 


To check our calculation for the self-energy we can re¬ 
late the RPT to the perturbation theory of Yamada and 
Yosida®] by replacing U e in Eq. (23) with the bare U, 


the parameters id, A with their bare counter-parts and 
setting all the counter-terms equal to zero. Similarly, 
we can check the four-vertex against the calculation of 
RefP. We find that the analytic results 


^E(O) 



P t 4.(0,0; 0,0) = U 




u 4 + ... 


(29) 


where u = U/ttA, are reproduced by our calculation. 


IV. RESULTS 


t/ = f n (0,0;0,0) = U e + J2 a nUe n - (24) 

n—2 

Working order-by-order we can invert this equation: 

00 

U e = U + Y J PnU n , (25) 

n—2 

where 

/?2 = — 0 - 2 , 

(3 3 = 2 a 9 — a 3 , 

/?4 = —5«2 + 5a 2 a3 — oq, 

/?5 = 14(3^2 — 21 G7<3 3 d“ 6(32*34 V 3(3g — (3 5 . (26) 


In this section we present numerical results for the ir¬ 
reducible self-energy and resultant spectral density. For 
all the calculations we fix 7 rA = H/100, where D is the 
conduction band width, and U = 3nA. In our discussion 
we consider the following parameter configurations: (i) a 
symmetric model with ed = —U/2; (ii) a model with 
weak asymmetry Ed = —1.27rA; and (iii) a model with 
pronounced asymmetry and = — 37 rA. In all cases, 
the scale of the problem is set by the renormalized den¬ 
sity of states at the Fermi level, 


Po 


A / 7T 

i d 2 + A 2 ' 


(30) 


We use po to define the Kondo temperature T K as 
Tk = 1 / 4 / 90 ; in the Kondo limit this reduces to the usual 
definition of Tk in terms of the susceptibilitjP We will 


Thus, the calculation of the self-energy in terms of U e 
can be rewritten as a series in U 

OO 

£(u >) = Y,6 n {u)U n , (27) 

n=2 


A. Self-energy 

For comparison purposes we will juxtapose the renor¬ 
malized self-energy obtained from RPT with the corre¬ 
sponding result obtained from the NRG. Since NRG cal¬ 
culations are set up with the bare self-energy in mind, we 
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have to use Eq. Q to relate the two quantities. For the 
imaginary part we find that 

ImS ff (w) = ^ImE^w). (31) 


We remark that inaccuracies in our NRG calculation re¬ 
sult in a slightly non-zero ImEo^O); to correct for this 
we offset our results by a small imaginary number. An 
equation similar to Eq. (311 can be derived for the real 


part; in practice however, it is of limited use. The reason 
is that the renormalized parameters are not determined 
from the NRG using the definitions in Eq. ©> but from 
the effective linear chain Hamiltonian (for more informa¬ 
tion see the Appendix of RefP^l). Since the low-energy 
limit of the real part of Eq. (|9j) relies crucially on the nu¬ 
merical cancellation between E(w) and (z— l)w, we found 
it difficult to obtain a result for ReE(w) with a vanishing 
derivative at zero frequency. It is for this reason that we 
only show NRG results for ImE(u). Note that due to the 
NRG’s successive elimination of higher energy scales we 
expect that its estimate for E(w) will only be accurate in 
the low frequency sector. 

We begin by discussing the symmetric model (case 
(i)) . As we have already remarked, in the presence of 
particle-hole symmetry id = 0 and the propagator is an 
odd function of frequency. The parity of the Green’s 
function also results in the cancellation of the particle- 
hole and particle-particle propagators 


W2 ) - _ ^.| 1 _ 1 j(wi,a;2), (32) 


and consequently the odd-order terms of the self-energy 
vanish for all frequencies. Using the method described in 
RefP 2 we find that A = 2.54 x 10 -4 and U = 7.95 x 10 -4 , 
giving U/ttA = 0.994. That this ratio is almost 1 is no 
coincidence; in the Kondo limit a universal scale, the 
Kondo temperature Tk emerges and U —>■ 4Tk,ttA — y 
4 Tj^. 

In all cases, the RPT estimate of the renormalized 
self-energy will always be exact in the limit uj —> 0, in 
the sense of Eq. ©>• At particle-hole symmetry, and for 
small but finite frequencies |wpo| 1, the w 2 coefficient 
of E(w) is reproduced exactlj^lby the second-order cal¬ 
culation; the fourth order term does not contribute terms 
of order uj 2 . As we increase |w| the dominant term is the 
w 4 term, which is not exactly given by the second-order 
calculation; this is corrected by the fourth-order contri¬ 
bution, extending the domain of validity of the RPT. At 
higher frequencies the disparity between the two RPT 
curves suggests the breakdown of the expansion; to ob¬ 
tain reliable results one must calculate the higher-order 
terms. We thus expect more elaborate approximations 
within RPT to continually extend the domain of validity 
of the resultant E(w). Were we to have, somehow, the 
ability to take into account all Feynman diagrams to all 
orders we would recover aE(w) exact for all frequencies, 
though clearly then we would also be able to calculate 
E(w) in the first place. 


Next, we discuss a slightly asymmetric model with = 
— 1.27rA, for which we find that id = 2.14 x 1CU 5 , A = 
2.93 x 10 -4 and U = 9.18 x 10 -4 . The non-zero but 
small in magnitude id now gives rise to odd-order terms 
in E(w). From Fig. ([4]) we see that the third and fifth 
order terms are small in value and essentially only slightly 
modify the second and fourth order curves respectively. 
To determine the stability of the series to order n we see 
that to simply look at the n + 1 term is not sufficient, 
since the n + 2 terms can still remain important. 

Finally, we turn our attention to Fig. ([5]) and case (iii), 
a very asymmetric model with = —37 tA. For small |w| 
we find again that the RPT is in good agreement with the 
NRG and that the lower order contributions dominate 
the result. A dramatic breakdown of the expansion, in 
both the real and imaginary components, is evident for 
small negative value of u> in contrast to cases (i) and (ii) 
where it is reasonably well-behaved even around wpo ~ 
0.5. 


B. Spectral densities 


We define the interacting Green’s function 

G a (u>) = -=-=—-, 

w — id + iA — E fl (w) 


(33) 


where Eis the retarded self-energy, which, as usual, 
differs from the causal E(w) one obtains from the dia- 
grammatics only by a sign(w) term in the imaginary part. 
We now turn our attention to the quasi-particle spectral 
density, defined as p(u>) = -ilmG(u) and express it in 
terms of the renormalized parameters as 


pH = - 


A — ImE fl (w) 


7r ( w - i d - ReE R {u>)) 2 + (A - ImE«(o;)) 2 

(34) 

Note that the spectral density is sensitive to the reducible 
self-energy. The low-frequency properties of the n’th or¬ 
der spectral density are thus the result of competition be¬ 
tween higher-order irreducible self-energies and powers of 
Ce 0 ! (w) combined with powers of lower-order irreducible 
self-energies. 

Results for cases (i), (ii) and (iii) are shown in Figs.^ 
and [8] respectively. To avoid the instability in the real 
part of E(w) we do not use Eq. (34) directly to extract 


the result from the NRG; instead we use the NRG result 


for p{uj) and, from Eq. (10), p(ui) = p(uj)/z. Note that 


Po is calculated from the renormalised parameters, and 
consequently, due to inaccuracies in our NRG calculation, 
the ratio p(ui = 0)/po m &y differ from 1. 

For the particle-hole symmetric model of case (i) we 
find again that the second and fourth order terms coin¬ 
cide for small uj and start deviating from each other and 
the NRG result as |w| is increased, with the RPT4 curve 
being in closer agreement with the latter. This picture 
persists in Fig. Q where, as in the case of the self-energy, 









FIG. 3: The real (left) and imaginary (right) parts of the renormalized self-energy for the particle-hole symmetric 

model. 




FIG. 4: The real (left) and imaginary (right) parts of the renormalized self-energy for the weakly asymmetric model 

(model (ii)). 


two pairs of similar curves emerge. We remark that at 
particle-hole symmetry p(w = 0) = (7 tA) - 1 ; away from 
particle-hole symmetry p(uj = 0) is given exactly by the 
renormalized parameters as per Eq. (30). 


V. CONCLUSION 

In this paper we have presented a relatively simple 
way of automating the calculation of the renormalized 
self-energy so that it can be carried out by a computer 
without user intervention. Our presentation was logically 
partitioned into three steps: the diagram generation, the 
application of the rules and analytic simplification of the 
integrals and the numerical integration of the diagrams. 

To illustrate the usefulness of the method we have cal¬ 
culated the self-energy up to fifth-order inclusive, a cal¬ 


culation which would otherwise be extremely tedious to 
perform by hand. We performed the calculation for three 
possible values of the asymmetry and found the results 
to be in good agreement with the NRG in the meaningful 
low frequency region. In all cases we find that the higher 
order term’s contribute more at higher frequencies, with 
RPT becoming increasingly more accurate as w —> 0. 

Though our discussion has been confined to the single¬ 
impurity Anderson model, the method here is readily 
generalisable to other models by replacing the Green’s 
function with the appropriate one. Unfortunately, the 
calculation of the fc-particle/hole propagator of the Ap¬ 
pendix relies on the linear dependence of [G^°l(u;)] _1 on 
w and thus does not generalize. This is not too serious 
an obstacle, for we can always numerically evaluate and 
tabulate the cases k = 2, 3 which most commonly appear. 
A further generalization can be realized by replacing the 
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FIG. 5: The real (left) and imaginary (right) parts of the renormalized self-energy for the strongly asymmetric 

model (model (iii)). 




FIG. 6: The spectral density for the particle-hole 
symmetric model (model (i)). 


FIG. 8: The spectral density for the very asymmetric 
model (model (iii)). 



scalar propagator with a matrix quantity to obtain calcu¬ 
late the behaviour of the impurity out of equilibrium 4 ^. 
This is of particular importance in the study of quantum 
dots under a bias voltage. 
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Appendix: Appendix 


In Eq. (141 we define a /c-particle/hole propagator as 

/ oo & 

d^n^w +wo. (at) 

-OO „•_i 


FIG. 7: The spectral density for model (ii), with some 
asymmetry. 
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where Si £ {—1,1}. Due to the sign term in the causal 
Green’s function (Eq. the integrand can be thought 
of as a piecewise function in a/. We begin by identifying 
the points j/i, which we will call nodes, where Siyi + tOi = 
0. We can ensure by appropriate labelling of the nodes 
that —oo > yi > 2/2 • > 2/n- For brevity, we write the 

i’th Green’s function in the form 


The decomposition of the real axis into intervals on which 
the integrand does not change form means that we can 
perform a partial fraction decomposition with coefficients 
specific to the region. Hence we can write 


1 

(u/ - ai)(w' - a 2 ) ■ ■ • (w' - a n ) 


E 


a 

u>' - O-i ’ 


(A.6) 


g£ ] M 


Si 

10 - ai(w) ’ 


(A.2) 


where ctj = Si [—oJi + ed i(Ti — iA sign(sjW + w*)]. We tem¬ 
porarily make the assumption, which will be lifted later, 
that all the a.i are distinct. Written as a product of terms 
of the form of Eq. (A.2), the integrand depends explicitly 
on to 1 11 * * but also implicitly through the dependence of the 
a. Having identified the nodes we can rewrite Eq. (A.l) 
as 


n— 1 

* _n n^(wi,...,w n ) = J(yi,y n )+^2F(yi,Vi+i), (A.3) 

2—1 


where 


F{yuyi+ 1 ) 


J(yi,y n ) 


fVi+l 

Jyi 


k 

d lo' G^] (siOj' + LOi) 

i=1 


[ ryi k 

lim < / du/ TT G ^ (siJ + w*) 

A ^oo [ J_ A f \ 



d to' G^j ( SiLO' + LOi) 


i—1 


(A.4) 

(A.5) 


where = 1 //'(%■), /(w) = (w - aq)... (w - a n ). We 
can now simply integrate each partial fraction separately. 
Let aW and / 3 W denote the values of the relevant quan¬ 
tities in the region (j/j, 2/i+i)- Then 


F(y„y l+1 )^±f>fLJ V ^p j 
i -i V t/i-o; / 

n 

J{yi,y n ) = E [^° )Ln _ a j 0) ) _ 

j=i 

o(n+l) T / (n+l)\ 

^ Ln ) 


(A.7) 


where Ln(z) denotes the principal branch of the complex 
logarithm defined as ln|z| + ?Arg(z), — n < Arg(z) < n 
and the labels 0, n+1 denote the values of the underlying 
quantities in the intervals (— 00 , 2 / 1 ) and (y n ,oo) respec¬ 
tively. 

In practical applications we may encounter numerical 
difficulties if the to, are such that any two ct; in a par¬ 
ticular region coincide, or nearly coincide This will cause 
our partial fraction to break down. We deal with this in 
a crude yet effective manner: when any on, aj are too 
close to each other, we separate them by a very small, 
arbitrary constant. After separating the offending Ui,a.j 
it is important to update the values of the corresponding 
2 /,;, y-j to ensure the consistency of the calculation. 
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